simulation results/results_pois_data/ratio of lambdas.R

In BriceonWiley/IntegratedLikelihood.R:

library(tidyverse); theme_set(theme_minimal())
library(parallel)
###############################################################################
calculate_mles <- function(data) {
  z1 <- data[1]; z2 <- data[2]; m01 <- data[3]; m02 <- data[4]
  y01 <- data[5]; y02 <- data[6]; N0 <- data[7]; N <- data[8]
  z01 <- m01 - y01 + y02; z02 <- m02 - y02 + y01
  lam1hat <- (z1 + m01 + (y01 * z2) / z02 - (y02 * z1) / z01) / (N + N0)
  lam2hat <- (z2 + m02 + (y02 * z1) / z01 - (y01 * z2) / z02) / (N + N0)
  theta1hat <- (y01 * z01 * (z2 + z02)) /
    (y01 * z01 * (z2 + z02) + (m01 - y01) * z02 * (z1 + z01))
  theta2hat <- (y02 * z02 * (z1 + z01)) /
    (y02 * z02 * (z1 + z01) + (m02 - y02) * z01 * (z2 + z02))
  c(
    phihat = lam1hat / lam2hat, 
    lam2hat = lam2hat,
    theta1hat = theta1hat,
    theta2hat = theta2hat # removed lam1hat
  )
}
Hessian <- function(params, data) {
  phi <- params[1]; lam2 <- params[2]; theta1 <- params[3]; theta2 <- params[4]
  z1 <- data[1]; z2 <- data[2]; m01 <- data[3]; m02 <- data[4]
  y01 <- data[5]; y02 <- data[6]; N0 <- data[7]; N <- data[8]
  matrix(
    c(
      -(m01 / phi^2) - 
        (z1 * (1 - theta1)^2) / (theta2 + (1 - theta1) * phi)^2 - 
        (z2 * theta1^2) / (1 - theta2 + theta1 * phi)^2, 
      -N - N0,
      (z1 * (1 - theta1) * phi) / (theta2 + (1 - theta1) * phi)^2 -
        z1 / (theta2 + (1 - theta1) * phi) - 
        (z2 * theta1 * phi) / (1 - theta2 + theta1 * phi)^2 + 
        z2 / (1 - theta2 + theta1 * phi), 
      -((z1 * (1 - theta1)) / (theta2 + (1 - theta1) * phi)^2) + 
        (z2 * theta1) / (1 - theta2 + theta1 * phi)^2, 
      -N - N0, 
      -((m01 + m02 + z1 + z2) / lam2^2), 
      0, 
      0, 
      (z1 * (1 - theta1) * phi) / (theta2 + (1 - theta1) * phi)^2 - 
        z1 / (theta2 + (1 - theta1) * phi) - 
        (z2 * theta1 * phi) / (1 - theta2 + theta1 * phi)^2 + 
        z2 / (1 - theta2 + theta1 * phi), 
      0, 
      -((m01 - y01) / (1 - theta1)^2) - 
        y01 / theta1^2 - 
        (z1 * phi^2) / (theta2 + (1 - theta1) * phi)^2 - 
        (z2 * phi^2) / (1 - theta2 + theta1 * phi)^2, 
      (z1 * phi) / (theta2 + (1 - theta1) * phi)^2 + 
        (z2 * phi) / (1 - theta2 + theta1 * phi)^2, 
      -((z1 * (1 - theta1)) / (theta2 + (1 - theta1) * phi)^2) + 
        (z2 * theta1) / (1 - theta2 + theta1 * phi)^2, 
      0, 
      (z1 * phi) / (theta2 + (1 - theta1) * phi)^2 + 
        (z2 * phi) / (1 - theta2 + theta1 * phi)^2, 
      -((m02 - y02) / (1 - theta2)^2) - 
        y02 / theta2^2 - 
        z1 / (theta2 + (1 - theta1) * phi)^2 - 
        z2 / (1 - theta2 + theta1 * phi)^2
    ), nrow = 4, ncol = 4
  )
}
Information <- function(params, data) {
  phi <- params[1]; lam2 <- params[2]; theta1 <- params[3]; theta2 <- params[4]
  z1 <- data[1]; z2 <- data[2]; m01 <- data[3]; m02 <- data[4]
  y01 <- data[5]; y02 <- data[6]; N0 <- data[7]; N <- data[8]
  matrix(
    c(
      # (1,1)
      N0 * lam2 / phi +
        ((1 - theta1)^2 * N * lam2) / (theta2 + (1 - theta1) * phi) +
        (theta1^2 * N * lam2) / (1 - theta2 + theta1 * phi),
      # (2,1)
      N + N0,
      # (3,1)
      -(1 - theta1) * phi * N * lam2 / (theta2 + (1 - theta1) * phi) +
        theta1 * N * phi * lam2 / (1 - theta2 + theta1 * phi),
      # (4,1)
      (1 - theta1) * N * lam2 / (theta2 + (1 - theta1) * phi) -
        theta1 * N * lam2 / (1 - theta2 + theta1 * phi),
      # (1,2)
      N + N0,
      # (2,2)
      (N0 + N) * (phi + 1) / lam2,
      # (3,2)
      0,
      # (4,2)
      0,
      # (1,3)
      -(1 - theta1) * phi * N * lam2 / (theta2 + (1 - theta1) * phi) +
        theta1 * N * phi * lam2 / (1 - theta2 + theta1 * phi),
      # (2,3)
      0,
      # (3,3)
      (N0 * phi * lam2 - m01 * theta1) / (1 - theta1)^2 +
        m01 / theta1 +
        phi^2 * N * lam2 / (theta2 + (1 - theta1) * phi) +
        phi^2 * N * lam2 / (1 - theta2 + theta1 * phi),
      # (4,3)
      - phi * N * lam2 / (theta2 + (1 - theta1) * phi) -
        phi * N * lam2 / (1 - theta2 + theta1 * phi),
      # (1,4)
      (1 - theta1) * N * lam2 / (theta2 + (1 - theta1) * phi) -
        theta1 * N * lam2 / (1 - theta2 + theta1 * phi),
      # (2,4)
      0,
      # (3,4)
      - phi * N * lam2 / (theta2 + (1 - theta1) * phi) -
        phi * N * lam2 / (1 - theta2 + theta1 * phi),
      # (4,4)
      (N0 * lam2 + m02 * theta2) / (1 - theta2)^2 +
        m02 / theta2 +
        N * lam2 / (theta2 + (1 - theta1) * phi) +
        N * lam2 / (1 - theta2 + theta1 * phi)
    ), nrow = 4, ncol = 4
  )
}
wald_phi_CI <- function(data, alpha = 0.05) {
  crit <- qnorm(1 - alpha / 2)
  mles <- calculate_mles(data)
  # if(any(mles == 0)) return(tibble(W_Lower = NaN, W_Upper = NaN))
  phihat <- mles[1]
  se <- tryCatch(
    solve(Information(mles, data))[1, 1],
    error = function(e) NaN
  )
  if (is.nan(se)) return(tibble(W_Lower = NaN, W_Upper = NaN))
  else {
    lower <- phihat - crit * sqrt(se)
    upper <- phihat + crit * sqrt(se)
    tibble(W_Lower = lower, W_Upper = upper)
  }
}
log_likelihood_phi <- function(params, phi, data) {
  lam2 <- params[1]; theta1 <- params[2]; theta2 <- params[3]
  z1 <- data[1]; z2 <- data[2]; m01 <- data[3]; m02 <- data[4]
  y01 <- data[5]; y02 <- data[6]; N0 <- data[7]; N <- data[8]
  
  m01 * log(phi * lam2) + m02 * log(lam2) - N0 * (phi * lam2 + lam2) +
    y01 * log(theta1) + (m01 - y01) * log(1 - theta1) +
    y02 * log(theta2) + (m02 - y02) * log(1 - theta2) +
    z1 * log(phi * lam2 * (1 - theta1) + lam2 * theta2) +
    z2 * log(lam2 * (1 - theta2) + phi * lam2 * theta1) - 
    N * (phi * lam2 + lam2)
}
score_function_phi <- function(params, phi, data) {
  lam2 <- params[1]; theta1 <- params[2]; theta2 <- params[3]
  z1 <- data[1]; z2 <- data[2]; m01 <- data[3]; m02 <- data[4]
  y01 <- data[5]; y02 <- data[6]; N0 <- data[7]; N <- data[8]
  
  -N * lam2 - N0 * lam2 + m01 / phi + 
    (z1 * (1 - theta1) * lam2) / (theta2 * lam2 + (1 - theta1) * lam2 * phi) + 
    (z2 * theta1 * lam2) / ((1 - theta2) * lam2 + theta1 * lam2 * phi)
}
score_vals_phi <- function(phi, data, inits, crit) {
  mles <- optim(
    inits[2:4], log_likelihood_phi, phi = phi, data = data, 
    lower = c(1.e-5, 1.e-5, 1.e-5), upper = c(100, 1 - 1.e-5, 1 - 1.e-5), 
    method = "L-BFGS-B", control = list(fnscale = -1)
  )$par
  val <- solve(Information(c(phi, mles), data))[1, 1]
  (score_function_phi(mles, phi, data))^2 * val - crit
}
score_phi_CI <- function(data, alpha = 0.05) {
  # browser()
  crit <- qchisq(1 - alpha, 1)
  inits <- calculate_mles(data)
  if(
    any(inits == 0) | any(is.nan(inits)) | 
    any(inits == Inf) | any(inits == -Inf)
  ) {
    return(tibble(S_Lower = NaN, S_Upper = NaN))
  }
  lower <- tryCatch(
    uniroot(
      score_vals_phi, c(1e-5, inits[1]), data = data,
      inits = inits, crit = crit
    )$root,
    error = function(e) NaN
  )
  upper <- tryCatch(
    uniroot(
      score_vals_phi, c(inits[1], 100), data = data, 
      inits = inits, crit = crit
    )$root,
    error = function(e) {
      tryCatch(
        uniroot(
          score_vals_phi, c(inits[1], 1000), data = data, 
          inits = inits, crit = crit
        )$root,
        error = function(e) NaN
      )
    }
  )
  tibble(S_Lower = lower, S_Upper = upper)
}
###############################################################################
likelihood_orig_phi <- function(params, phi, data) {
  lam2 <- params[1, ]; theta1 <- params[2, ]; theta2 <- params[3, ]
  z1 <- data[1]; z2 <- data[2]; m01 <- data[3]; m02 <- data[4]
  y01 <- data[5]; y02 <- data[6]; N0 <- data[7]; N <- data[8]
  
  log_lik <- 
    m01 * log(phi) + (m01 + m02 + z1 + z2) * log(lam2) - 
    lam2 * (phi + 1) * (N0 + N) +
    z1 * log(phi * (1 - theta1) + theta2) +
    z2 * log((1 - theta2) + phi * theta1) +
    y01 * log(theta1) + (m01 - y01) * log(1 - theta1) +
    y02 * log(theta2) + (m02 - y02) * log(1 - theta2) +
    0.001 * log(0.001) - lgamma(0.001) +
    (0.001 - 1) * log(lam2) - 0.001 * lam2

  matrix(exp(log_lik), ncol = ncol(params))
}
likelihood_phi <- function(phi, data) {
  hcubature(
    likelihood_orig_phi, lowerLimit = c(0, 0, 0), upperLimit = c(Inf, 1, 1),
    phi = phi, data = data, vectorInterface = TRUE
  )$integral
}
closed_phi <- function(phi, data) {
  # browser()
  z1 <- data[1]; z2 <- data[2]; m01 <- data[3]; m02 <- data[4]
  y01 <- data[5]; y02 <- data[6]; N0 <- data[7]; N <- data[8]
  
  first <- -(m01 + m02 + z1 + z2 + 0.001) *
    log((phi + 1) * (N0 + N) + 0.001)
  indices <- expand.grid(i = 0:z1, j = 0:z2)
  i <- indices$i; j <- indices$j
  combs <- choose(z1, i) * choose(z2, j)
  delt <- phi^(m01 + i + z2 - j)
  betas <- beta(m01 - y01 + i + 1, y01 + z2 - j + 1) *
    beta(m02 - y02 + j + 1,  y02 + z1 - i + 1)
  second <- combs * delt * betas
  exp(first) * sum(second)
}
likelihood_ratio_phi <- function(phi, H_1, data, crit, closed) {
  if (closed) H_0 <- closed_phi(phi, data)
  else H_0 <- likelihood_phi(phi, data)
  -2 * log(H_0 / H_1) - crit
}
intLik_phi_CI <- function(data, alpha = 0.05, closed_form = TRUE) {
  # browser()
  crit <- qchisq(1 - alpha, 1)
  if(closed) opt <- optimize(closed_phi, c(1e-6, 10), data = data, maximum = T)
  else opt <- optimize(likelihood_phi, c(1e-6, 10), data = data, maximum = T)
  phihat <- opt$maximum; H_1 <- opt$objective
  
  # phihat <- calculate_mles(data)[1]
  # if (closed) H_1 <- closed_phi(phihat, data)
  # else H_1 <- likelihood_phi(phihat, data)
  lower <- tryCatch(
    uniroot(
      likelihood_ratio_phi, c(1e-6, phihat),
      H_1 = H_1, data = data, crit = crit, closed = closed
    )$root,
    error = function(e) NaN
  )
  upper <- tryCatch(
    uniroot(
      likelihood_ratio_phi, c(phihat, 10),
      H_1 = H_1,
      data = data, crit = crit, closed = closed
    )$root,
    error = function(e) {
      tryCatch(
        uniroot(
          likelihood_ratio_phi, c(phihat, 100),
          H_1 = H_1,
          data = data, crit = crit, closed = closed
        )$root,
        error = function(e) {
          tryCatch(
            uniroot(
              likelihood_ratio_phi, c(phihat, 1000),
              H_1 = H_1,
              data = data, crit = crit, closed = closed
            )$root,
            error = function(e) NaN
          )
        }
      )
    }
  )
  tibble(ILR_Lower = lower, ILR_Upper = upper)
}
simulate_phi <- function(N0, N, theta1, theta2 = theta1, lam1 = 20, lam2 = 15, 
                         len = 10000, alpha = 0.05, adjustment = 0,
                         closed, summary = FALSE) {
    # browser()
    samples <- rerun(len, {
      m01 <- rpois(1, N0 * lam1)
      m02 <- rpois(1, N0 * lam2)
      y01 <- rbinom(1, m01, theta1)
      y02 <- rbinom(1, m02, theta2)
      z1 <- rpois(1, N * (lam1 * (1 - theta1) + lam2 * theta2))
      z2 <- rpois(1, N * (lam2 * (1 - theta2) + lam1 * theta1))
      if (m01 == 0) m01 <- adjustment
      if (m02 == 0) m02 <- adjustment
      if (y01 == 0) y01 <- adjustment
      if (y02 == 0) y02 <- adjustment
      tibble(z1 = z1, z2 = z2, m01 = m01, m02 = m02, y01 = y01, y02 = y02)
    }) %>% 
      bind_rows()
    tryCatch(
      walds <- pmap_dfr(
        samples, ~ wald_phi_CI(c(..1, ..2, ..3, ..4, ..5, ..6, N0, N), alpha)
      ),
      error = function(e) browser()
    )
    scores <- pmap_dfr(
      samples, ~ score_phi_CI(c(..1, ..2, ..3, ..4, ..5, ..6, N0, N), alpha)
    )
    ILR <- pmap_dfr(
      samples, ~ intLik_phi_CI(
        c(..1, ..2, ..3, ..4, ..5, ..6, N0, N), alpha, closed
      )
    )
    intervals <- bind_cols(walds, scores, ILR)
    NaNs <- intervals %>% 
      pmap_lgl(~ any(is.nan(c(..1, ..2, ..3, ..4, ..5, ..6)))) %>% 
      sum()
    values <- c(
      "lambda1" = lam1, "lambda2" = lam2, "phi" = round(lam1 / lam2, 4),
      "theta1" = theta1, "theta2" = theta2, "N0" = N0, "N" = N, 
      "B" = len, "NaN" = NaNs
    )
    results <- list(
      "values" = values,
      "samples" = samples,
      "intervals" = intervals
    )
    
    # producing final product
    if (summary) {
      coverages <- coverage_phi(results, ignore_NaN = TRUE)
      widths <- width_phi(results, ignore_NaN = TRUE)
      summarized <-
        bind_rows(coverages, widths) %>%
        select(-N0, -N, -phi, -theta1, -theta2) %>%
        t() %>%
        `colnames<-`(c("Coverage Probability", "Average Width"))
      message(
        glue("Removed {results$values['NaN']} bad intervals in summary calculations.")
      )
      list(
        "values" = results$values,
        "summary" = round(summarized, 4)
      )
    } else results
}
coverage_phi <- function(x, ignore_NaN = FALSE) {
  N0 <- x$values["N0"]; N <- x$values["N"]
  phi <- x$values["phi"]; 
  theta1 <- x$values["theta1"]; theta2 <- x$values["theta2"]
  int_names <- names(x$intervals)
  if(ignore_NaN) {
    usable <- filter(x$intervals, !is.nan(W_Lower)) %>% 
      filter(!is.nan(S_Upper)) %>% 
      filter(!is.nan(ILR_Upper))
    lower <- select(usable, str_subset(int_names, "Lower"))
    upper <- select(usable, str_subset(int_names, "Upper"))
  } else {
    lower <- select(x$intervals, str_subset(int_names, "Lower"))
    upper <- select(x$intervals, str_subset(int_names, "Upper"))
  }
  coverages <- map2(lower, upper, ~ map2(.x, .y, ~ between(phi, .x, .y))) %>%
    transpose() %>%
    bind_rows() %>%
    summarise_all(mean)
  names(coverages) <- str_remove(names(coverages), "[_]?Lower[_]?")
  coverages %>%
    add_column(
      "N0" = N0, "N" = N, "phi" = phi, 
      "theta1" = theta1, "theta2" = theta2, .before = TRUE
    )
}
width_phi <- function(x, ignore_NaN = FALSE) {
  N0 <- x$values["N0"]; N <- x$values["N"]
  phi <- x$values["phi"]; 
  theta1 <- x$values["theta1"]; theta2 <- x$values["theta2"]
  int_names <- names(x$intervals)
  if (ignore_NaN) {
    usable <- filter(x$intervals, !is.nan(W_Lower)) %>% 
      filter(!is.nan(S_Upper)) %>% 
      filter(!is.nan(ILR_Upper))
    lower <- select(usable, str_subset(int_names, "Lower"))
    upper <- select(usable, str_subset(int_names, "Upper"))
  } else {
    lower <- select(x$intervals, str_subset(int_names, "Lower"))
    upper <- select(x$intervals, str_subset(int_names, "Upper"))
  }
  widths <- map2(lower, upper, ~ map2(.x, .y, ~ (.y - .x))) %>%
    transpose() %>%
    bind_rows() %>%
    summarise_all(mean)
  names(widths) <- str_remove(names(widths), "[_]?Lower[_]?")
  widths %>%
    add_column(
      "N0" = N0, "N" = N, "phi" = phi, 
      "theta1" = theta1, "theta2" = theta2, .before = TRUE
    )
}
###############################################################################
mix <- expand_grid(lam1 = 1:8, theta = seq(0.05, 0.95, 0.05))

phi_N0_1_N_2 <- mcMap(
  function(.l, .t) {
    simulate_phi(
      N0 = 1, N = 2, theta1 = .t, lam1 = .l, lam2 = 4, closed = T
    )
  }, .l = mix$lam1, .t = mix$theta, mc.cores = 4
)
phi_N0_1_N_3 <- mcMap(
  function(.l, .t) {
    simulate_phi(
      N0 = 1, N = 3, theta1 = .t, lam1 = .l, lam2 = 4, closed = T
    )
  }, .l = mix$lam1, .t = mix$theta, mc.cores = 4
)
phi_N0_2_N_4 <- mcMap(
  function(.l, .t) {
    simulate_phi(
      N0 = 2, N = 4, theta1 = .t, lam1 = .l, lam2 = 4, closed = T
    )
  }, .l = mix$lam1, .t = mix$theta, mc.cores = 4
)
phi_N0_2_N_6 <- mcMap(
  function(.l, .t) {
    simulate_phi(
      N0 = 2, N = 6, theta1 = .t, lam1 = .l, lam2 = 4, closed = T
    )
  }, .l = mix$lam1, .t = mix$theta, mc.cores = 4
)
phi_N0_3_N_6 <- mcMap(
  function(.l, .t) {
    simulate_phi(
      N0 = 3, N = 6, theta1 = .t, lam1 = .l, lam2 = 4, closed = T
    )
  }, .l = mix$lam1, .t = mix$theta, mc.cores = 4
)
phi_N0_3_N_9 <- mcMap(
  function(.l, .t) {
    simulate_phi(
      N0 = 3, N = 9, theta1 = .t, lam1 = .l, lam2 = 4, closed = T
    )
  }, .l = mix$lam1, .t = mix$theta, mc.cores = 4
)
phi_N0_4_N_8 <- mcMap(
  function(.l, .t) {
    simulate_phi(
      N0 = 4, N = 8, theta1 = .t, lam1 = .l, lam2 = 4, closed = T
    )
  }, .l = mix$lam1, .t = mix$theta, mc.cores = 4
)

###############################################################################
# wald_CI <- function(data, alpha = 0.05) {
#   crit <- qnorm(1 - alpha / 2)
#   mles <- unlist(calculate_mles(data))
#   lam1hat <- mles[1]
#   se <- solve(Information(mles, data))[1, 1]
#   lower <- lam1hat - crit * sqrt(se)
#   upper <- lam1hat + crit * sqrt(se)
#   tibble(W_Lower = lower, W_Upper = upper)
# }
# log_likelihood <- function(params, lam1, data) {
#   lam2 <- params[1]; theta1 <- params[2]; theta2 <- params[3]
#   z1 <- data[1]; z2 <- data[2]; m01 <- data[3]; m02 <- data[4]
#   y01 <- data[5]; y02 <- data[6]; N0 <- data[7]; N <- data[8]
#   
#   m01 * log(lam1) + m02 * log(lam2) - N0 * (lam1 + lam2) +
#     y01 * log(theta1) + (m01 - y01) * log(1 - theta1) +
#     y02 * log(theta2) + (m02 - y02) * log(1 - theta2) +
#     z1 * log(lam1 * (1 - theta1) + lam2 * theta2) +
#     z2 * log(lam2 * (1 - theta2) + lam1 * theta1) - N * (lam1 + lam2)
# }
# score_function <- function(params, lam1, data) {
#   lam2 <- params[1]; theta1 <- params[2]; theta2 <- params[3]
#   z1 <- data[1]; z2 <- data[2]; m01 <- data[3]; m02 <- data[4]
#   y01 <- data[5]; y02 <- data[6]; N0 <- data[7]; N <- data[8]
#   
#   -N - N0 + m01 / lam1 + 
#     (z2 * theta1) / (theta1 * lam1 + (1 - theta2) * lam2) + 
#     (z1 * (1 - theta1)) / ((1 - theta1) * lam1 + theta2 * lam2)
# }
# score_vals <- function(lam1, data, inits, crit) {
#   mles <- optim(
#     inits[-1], log_likelihood, lam1 = lam1, 
#     data = data, control = list(fnscale = -1)
#   )$par
#   val <- solve(Information(c(lam1, mles), data))[1, 1]
#   (score_function(mles, lam1, data))^2 * val - crit
# }
# score_CI <- function(data, alpha = 0.05) {
#   crit <- qchisq(1 - alpha, 1)
#   inits <- unlist(calculate_mles(data))
#   lower <- uniroot(
#     score_vals, c(1e-5, inits[1]), data = data,
#     inits = inits, crit = crit
#   )$root
#   upper <- uniroot(
#     score_vals, c(inits[1], 1000), data = data, 
#     inits = inits, crit = crit
#   )$root
#   tibble(S_Lower = lower, S_Upper = upper)
# }
# profile_ratio <- function(lam1, inits, H_1, data, crit) {
#   pars <- optim(
#     inits, log_likelihood, lam1 = lam1, 
#     data = data, control = list(fnscale = -1)
#   )$par
#   H_0 <- log_likelihood(pars, lam1, data)
#   -2 * (H_0 - H_1) - crit
# }
# profile_CI <- function(data, alpha = 0.05) {
#   # browser()
#   crit <- qchisq(1 - alpha, 1)
#   mles <- unlist(calculate_mles(data))
#   H_1 <- log_likelihood(mles[-1], mles[1], data)
#   lower <- uniroot(
#     profile_ratio, c(1e-5, mles[1]), inits = mles[-1],
#     H_1 = H_1, data = data, crit = crit
#   )$root
#   upper <- uniroot(
#     profile_ratio, c(mles[1], 100), inits = mles[-1],
#     H_1 = H_1, data = data, crit = crit
#   )$root
#   tibble(P_Lower = lower, L_Upper = upper)
# }
# profile_ratio_phi <- function(phi, inits, H_1, data, crit) {
#   # pars <- optim(
#   #   inits, log_likelihood_phi, phi = phi,
#   #   data = data, control = list(fnscale = -1)
#   # )$par
#   pars <- optim(
#     inits, log_likelihood_phi, phi = phi, data = data, 
#     lower = c(1.e-5, 1.e-5, 1.e-5), upper = c(100, 1 - 1.e-5, 1 - 1.e-5), 
#     method = "L-BFGS-B", control = list(fnscale = -1)
#   )$par
#   H_0 <- log_likelihood_phi(pars, phi, data)
#   -2 * (H_0 - H_1) - crit
# }
# profile_phi_CI <- function(data, alpha = 0.05) {
#   # browser()
#   crit <- qchisq(1 - alpha, 1)
#   mles <- calculate_mles(data)
#   if(any(mles == 0)) return(tibble(P_Lower = NaN, P_Upper = NaN))
#   H_1 <- log_likelihood_phi(mles[2:4], mles[5], data)
#   lower <- tryCatch(
#     uniroot(
#       profile_ratio_phi, c(1e-5, mles[5]), inits = mles[2:4],
#       H_1 = H_1, data = data, crit = crit
#     )$root,
#     error = function(e) NaN
#   )
#   upper <- tryCatch(
#     uniroot(
#       profile_ratio_phi, c(mles[5], 100), inits = mles[2:4],
#       H_1 = H_1, data = data, crit = crit
#     )$root,
#     error = function(e) NaN
#   )
#   tibble(P_Lower = lower, P_Upper = upper)
# }
# Hessian <- function(params, data) {
#   lam1 <- params[1]; lam2 <- params[2]; theta1 <- params[3]; theta2 <- params[4]
#   z1 <- data[1]; z2 <- data[2]; m01 <- data[3]; m02 <- data[4]
#   y01 <- data[5]; y02 <- data[6]; N0 <- data[7]; N <- data[8]
#   matrix(
#     c(
#       -(m01 / lam1^2) - 
#         (z2 * theta1^2) / (theta1 * lam1 + (1 - theta2) * lam2)^2 - 
#         (z1  * (1 - theta1)^2) / ((1 - theta1) * lam1 + theta2 * lam2)^2,
#       -((z2 * theta1 * (1 - theta2)) / (theta1 * lam1 + (1 - theta2) * lam2)^2) -
#         (z1 * (1 - theta1) * theta2) / ((1 - theta1) * lam1 + theta2 * lam2)^2,
#       -((z2 * theta1 * lam1) / (theta1 * lam1 + (1 - theta2) * lam2)^2) + 
#         z2 / (theta1 * lam1 + (1 - theta2) * lam2) + 
#         (z1 * (1 - theta1) * lam1) / ((1 - theta1) * lam1 + theta2 * lam2)^2 - 
#         z1 / ((1 - theta1) * lam1 + theta2 * lam2),
#       (z2 * theta1 * lam2) / (theta1 * lam1 + (1 - theta2) * lam2)^2 - 
#         (z1 * (1 - theta1) * lam2) / ((1 - theta1) * lam1 + theta2 * lam2)^2,
#       -((z2 * theta1 * (1 - theta2)) / (theta1 * lam1 + (1 - theta2) * lam2)^2) -
#         (z1 * (1 - theta1) * theta2) / ((1 - theta1) * lam1 + theta2 * lam2)^2,
#       -(m02 / lam2^2) - 
#         (z2 * (1 - theta2)^2) / (theta1 * lam1 + (1 - theta2) * lam2)^2 - 
#         (z1 * theta2^2) / ((1 - theta1) * lam1 + theta2 * lam2)^2, 
#       -((z2 * (1 - theta2) * lam1) / (theta1 * lam1 + (1 - theta2) * lam2)^2) + 
#         (z1 * theta2 * lam1) / ((1 - theta1) * lam1 + theta2 * lam2)^2, 
#       (z2 * (1 - theta2) * lam2) / (theta1 * lam1 + (1 - theta2) * lam2)^2 - 
#         z2 / (theta1 * lam1 + (1 - theta2) * lam2) - 
#         (z1 * theta2 * lam2) / ((1 - theta1) * lam1 + theta2 * lam2)^2 + 
#         z1/((1 - theta1) * lam1 + theta2 * lam2),
#       -((z2 * theta1 * lam1) / (theta1 * lam1 + (1 - theta2) * lam2)^2) + 
#         z2 / (theta1 * lam1 + (1 - theta2) * lam2) + 
#         (z1 * (1 - theta1) * lam1) / ((1 - theta1) * lam1 + theta2 * lam2)^2 - 
#         z1 / ((1 - theta1) * lam1 + theta2 * lam2), 
#       -((z2 * (1 - theta2) * lam1) / (theta1 * lam1 + (1 - theta2) * lam2)^2) + 
#         (z1 * theta2 * lam1) / ((1 - theta1) * lam1 + theta2 * lam2)^2, 
#       -((m01 - y01) / (1 - theta1)^2) - 
#         y01 / theta1^2 - 
#         (z2 * lam1^2) / (theta1 * lam1 + (1 - theta2) * lam2)^2 - 
#         (z1 * lam1^2) / ((1 - theta1) * lam1 + theta2 * lam2)^2, 
#       (z2 * lam1 * lam2) / (theta1 * lam1 + (1 - theta2) * lam2)^2 + 
#         (z1 * lam1 * lam2) / ((1 - theta1) * lam1 + theta2 * lam2)^2,
#       (z2 * theta1 * lam2) / (theta1 * lam1 + (1 - theta2) * lam2)^2 - 
#         (z1 * (1 - theta1) * lam2) / ((1 - theta1) * lam1 + theta2 * lam2)^2,
#       (z2 * (1 - theta2) * lam2) / (theta1 * lam1 + (1 - theta2) * lam2)^2 - 
#         z2 / (theta1 * lam1 + (1 - theta2) * lam2) -
#         (z1 * theta2 * lam2) / ((1 - theta1) * lam1 + theta2 * lam2)^2 + 
#         z1 / ((1 - theta1) * lam1 + theta2 * lam2), 
#       (z2 * lam1 * lam2) / (theta1 * lam1 + (1 - theta2) * lam2)^2 + 
#         (z1 * lam1 * lam2) / ((1 - theta1) * lam1 + theta2 * lam2)^2, 
#       -((m02 - y02) / (1 - theta2)^2) - 
#         y02 / theta2^2 - 
#         (z2 * lam2^2) / (theta1 * lam1 + (1 - theta2) * lam2)^2 - 
#         (z1 * lam2^2) / ((1 - theta1) * lam1 + theta2 * lam2)^2
#     ), nrow = 4, ncol = 4
#   )
# }
# Information <- function(params, data) {
#   lam1 <- params[1]; lam2 <- params[2]; theta1 <- params[3]; theta2 <- params[4]
#   z1 <- data[1]; z2 <- data[2]; m01 <- data[3]; m02 <- data[4]
#   y01 <- data[5]; y02 <- data[6]; N0 <- data[7]; N <- data[8]
#   matrix(
#     c(
#       # (1,1)
#       N0 / lam1 + 
#         (1 - theta1)^2 * N / ((1 - theta1) * lam1 + theta2 * lam2) +
#         theta1^2 * N / (theta1 * lam1 + (1 - theta2) * lam2),
#       # (1,2)
#       (1 - theta1) * theta2 * N / ((1 - theta1) * lam1 + theta2 * lam2) +
#         theta1 * (1 - theta2) * N / (theta1 * lam1 + (1 - theta2) * lam2),
#       # (1,3)
#       -(1 - theta1) * lam1 * N / ((1 - theta1) * lam1 + theta2 * lam2) +
#         theta1 * lam1 * N / (theta1 * lam1 + (1 - theta2) * lam2),
#       # (1,4)
#       -theta1 * lam2 * N / (theta1 * lam1 + (1 - theta2) * lam2) +
#         (1 - theta1) * lam2 * N / ((1 - theta1) * lam1 + theta2 * lam2),
#       # (2,1)
#       (1 - theta1) * theta2 * N / ((1 - theta1) * lam1 + theta2 * lam2) +
#         theta1 * (1 - theta2) * N / (theta1 * lam1 + (1 - theta2) * lam2),
#       # (2,2)
#       N0 / lam2 + 
#         theta2^2 * N / ((1 - theta1) * lam1 + theta2 * lam2) +
#         (1 - theta2)^2 * N / (theta1 * lam1 + (1 - theta2) * lam2),
#       # (2,3)
#       -theta2 * lam1 * N / ((1 - theta1) * lam1 + theta2 * lam2) +
#         (1 - theta2) * lam1 * N / (theta1 * lam1 + (1 - theta2) * lam2),
#       # (2,4)
#       theta2 * lam2 * N / ((1 - theta1) * lam1 + theta2 * lam2) -
#         (1 - theta2) * lam2 * N / (theta1 * lam1 + (1 - theta2) * lam2),
#       # (3,1)
#       -(1 - theta1) * lam1 * N / ((1 - theta1) * lam1 + theta2 * lam2) +
#         theta1 * lam1 * N / (theta1 * lam1 + (1 - theta2) * lam2),
#       # (3,2)
#       -theta2 * lam1 * N / ((1 - theta1) * lam1 + theta2 * lam2) +
#         (1 - theta2) * lam1 * N / (theta1 * lam1 + (1 - theta2) * lam2),
#       # (3,3)
#       N0 * lam1 / (1 - theta1)^2 - 
#         m01 * theta1 * (theta1^2 - (1 - theta1)^2) / ((1 - theta1)^2 * theta1^2) + 
#         lam1^2 * N / ((1 - theta1) * lam1 + theta2 * lam2) +
#         lam1^2 * N / (theta1 * lam1 + (1 - theta2) * lam2),
#       # (3,4)
#       -lam1 * lam2 * N / ((1 - theta1) * lam1 + theta2 * lam2) -
#         lam1 * lam2 * N / (theta1 * lam1 + (1 - theta2) * lam2),
#       # (4,1)
#       -theta1 * lam2 * N / (theta1 * lam1 + (1 - theta2) * lam2) +
#         (1 - theta1) * lam2 * N / ((1 - theta1) * lam1 + theta2 * lam2),
#       # (4,2)
#       theta2 * lam2 * N / ((1 - theta1) * lam1 + theta2 * lam2) -
#         (1 - theta2) * lam2 * N / (theta1 * lam1 + (1 - theta2) * lam2),
#       # (4,3)
#       -lam1 * lam2 * N / ((1 - theta1) * lam1 + theta2 * lam2) -
#         lam1 * lam2 * N / (theta1 * lam1 + (1 - theta2) * lam2),
#       # (4,4)
#       N0 * lam2 / (1 - theta2)^2 +
#         m02 * theta2 * (theta2^2 - (1 - theta2)^2) / ((1 - theta2)^2 * theta2^2) +
#         lam2^2 * N / ((1 - theta1) * lam1 + theta2 * lam2) +
#         lam2^2 * N / (theta1 * lam1 + (1 - theta2) * lam2)
#     ), nrow = 4, ncol = 4
#   )
# }
# Jacobian <- function(phi, lam2) {
#   J <- diag(c(lam2, 1, 1, 1))
#   J[1, 2] <- phi
#   J
# }
# plot_coverage_theta <- function(x, val) {
#   x %>% 
#     select(-N0, -N) %>% 
#     gather("Interval", "Coverage", -phi, -adj_val) %>% 
#     ggplot(aes(phi, Coverage)) +
#     geom_line(aes(color = Interval)) +
#     geom_hline(yintercept = 0.95, linetype = "dashed") +
#     facet_wrap(~ adj_val) +
#     ggtitle(glue("Theta = {val}"))
# }
# plot_width_theta <- function(x, val) {
#   x %>% 
#     select(-N0, -N) %>% 
#     gather("Interval", "Width", -phi, -adj_val) %>% 
#     ggplot(aes(phi, Width)) +
#     geom_line(aes(color = Interval)) +
#     facet_wrap(~ adj_val) +
#     ggtitle(glue("Theta = {val}"))
# }
# phi_N0_1_N_2_adj <- mcMap(
#   function(.l, .t) {
#     simulate_phi(
#       N0 = 1, N = 2, theta1 = .t, lam1 = .l, lam2 = 4, len = 2500, closed = T,
#       adjustment = 0.5
#     )
#   }, .l = mix$lam1, .t = mix$theta, mc.cores = 4
# )
# phi_N0_1_N_3_adj <- mcMap(
#   function(.l, .t) {
#     simulate_phi(
#       N0 = 1, N = 3, theta1 = .t, lam1 = .l, lam2 = 4, len = 2500, closed = T,
#       adjustment = 0.5
#     )
#   }, .l = mix$lam1, .t = mix$theta, mc.cores = 4
# )
# phi_N0_2_N_4_adj <- mcMap(
#   function(.l, .t) {
#     simulate_phi(
#       N0 = 2, N = 4, theta1 = .t, lam1 = .l, lam2 = 4, len = 2500, closed = T,
#       adjustment = 0.5
#     )
#   }, .l = mix$lam1, .t = mix$theta, mc.cores = 4
# )
# phi_N0_2_N_6_adj <- mcMap(
#   function(.l, .t) {
#     simulate_phi(
#       N0 = 2, N = 6, theta1 = .t, lam1 = .l, lam2 = 4, len = 2500, closed = T,
#       adjustment = 0.5
#     )
#   }, .l = mix$lam1, .t = mix$theta, mc.cores = 4
# )
# phi_N0_3_N_6_adj <- mcMap(
#   function(.l, .t) {
#     simulate_phi(
#       N0 = 3, N = 6, theta1 = .t, lam1 = .l, lam2 = 4, len = 2500, closed = T,
#       adjustment = 0.5
#     )
#   }, .l = mix$lam1, .t = mix$theta, mc.cores = 4
# )
# phi_N0_3_N_9_adj <- mcMap(
#   function(.l, .t) {
#     simulate_phi(
#       N0 = 3, N = 9, theta1 = .t, lam1 = .l, lam2 = 4, len = 2500, closed = T,
#       adjustment = 0.5
#     )
#   }, .l = mix$lam1, .t = mix$theta, mc.cores = 4
# )
# phi_N0_4_N_8_adj <- mcMap(
#   function(.l, .t) {
#     simulate_phi(
#       N0 = 4, N = 8, theta1 = .t, lam1 = .l, lam2 = 4, len = 2500, closed = T,
#       adjustment = 0.5
#     )
#   }, .l = mix$lam1, .t = mix$theta, mc.cores = 4
# )